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Abstract. In this paper, we investigate the fundamental limits on how the inter- 
spike time of a neuron oscillator can be perturbed by the application of a bounded 
external control input (a current stimulus) with zero net electric charge accumulation. 
We use phase models to study the dynamics of neurons and derive charge-balanced 

■ controls that achieve the minimum and maximum inter-spike times for a given bound 
C$ • on the control amplitude. Our derivation is valid for any arbitrary shape of the phase 

response curve and for any value of the given control amplitude bound. In addition, we 
characterize the change in the structures of the charge-balanced time-optimal controls 
with the allowable control amplitude. We demonstrate the applicability of the derived 
^ 1 optimal control laws by applying them to mathematically ideal and experimentally 

■ observed neuron phase models, including the widely-studied Hodgkin-Huxley phase 

model, and by verifying them with the corresponding original full state-space models. 

■^j- ■ This work addresses a fundamental problem in the field of neural control and provides 

a theoretical investigation to the optimal control of oscillatory systems. 
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1. Introduction 

Neurons exhibit short-lasting voltage spikes known as action potentials, which are 
sensitive to external current stimuli [1]. The inter-spike time interval of a neuron 
characterizes its properties and can be controlled by use of external stimuli. The ability 
to control neuron spiking activities is fundamental to theoretical neuroscience, and the 
concept of effective control of such neurological behavior has led to the development of 
innovative therapeutic procedures [2, 3] for neurological disorders including deep brain 
stimulation (DBS) for Parkinson's disease and essential tremor [4, 5], where electrical 
pulses are used to inhibit pathological synchrony among neuron populations. In such 
neurological treatments and other applications such as the design of artificial cardiac 
pacemakers [6], it is of clinical importance to avoid long and strong electrical pulses in 
order to prevent the tissue from damage, as well as to maintain zero net electric charge 
accumulation over each stimulation cycle in order to suppress undesirable side effects. 
High levels of electric charge accumulation may trigger irreversible electrochemical 
reactions resulting in damage of neural tissues and corrosion of electrodes [7]. 

Motivated by these practical needs, in this paper we study the design of time- 
optimal controls for spiking neurons, which lead to the minimum and maximum inter- 
spike times and remain charge-balanced. We study the dynamics of neuron oscillators 
through phase models which are simplified yet accurate models that capture essential 
overall properties of an oscillating neuron [1, 8], and which form a standard nonlinear 
control system that characterizes the evolution of an oscillating system by a single 
variable, namely, the phase. Phase models are conventionally used to investigate the 
synchronization patterns and study the dynamical responses of oscillators where the 
inputs to the oscillatory systems are initially defined [8, 9, 10]. Recently, control- 
theoretic approaches, including calculus of variations and the maximum principle, have 
been employed to design external stimuli for optimal manipulation of the dynamic 
behavior of neuron oscillators. These include the design of minimum-power controls 
for spiking a single neuron at specified time instances [11, 12, 13], optimal waveforms 
for entrainment of neuron ensembles [14, 15, 16], and open- loop controls for establishing 
and maintaining a desired phase configuration, such as anti-phase for two coupled neuron 
oscillators [17]. Work on considering stochastic effects to neuron systems such as the 
optimal control of neuronal spiking activity receiving a class of random synaptic inputs 
has also been investigated [18]. In addition, controllability of an ensemble of uncoupled 
neurons was explored for various mathematically ideal phase models, where an effective 
computational optimal control method based on pseudospectral approximations was 
employed to construct optimal controls that elicit simultaneous spikes of a neuron 
ensemble [19, 20]. The derivation of time-optimal and spike timing controls for spiking 
neurons has been attempted for limited classes of control functions [21, 22], however, 
a complete characterization of the optimal solutions has not been provided, and an 
analytical and systematic approach for synthesizing the time-optimal controls has been 
missing. 
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In this paper, we derive charge-balanced time-optimal controls for a given bound on 
the control amplitude and fully characterize the possible range of neuron spiking times 
determined by such optimal controls. Employing techniques from the optimal control 
theory, we are able to reveal different structures of the time-optimal controls that vary 
with the allowable bound of the control amplitude. Moreover, we validate these controls 
derived according to phase models by applying them to the corresponding original full 
state-space neuron models. As a demonstration, the validation is performed using the 
Hodgkin- Huxley equations [23], where the spiking behavior of the state-space model 
shows great qualitative agreement with that of the phase model and which demonstrates 
the applicability of our theoretical results based on the phase model. Such an important 
validation, which is largely lacking in the literature, allows us to explore the fundamental 
limits of the phase model as an approximation of state-space models. 

This paper is organized as follows. In Section 2, we consider the time-optimal 
control of a general phase oscillator and derive the charge-balanced minimum-time and 
maximum-time controls with constrained control amplitude by using the Pontryagin's 
maximum principle [24]. In Section 3, we apply the derived optimal control strategies 
to both mathematically ideal and experimentally observed phase models, including the 
well-known SNIPER [8], Hodgkin-Huxley, and Morris-Lecar [25] phase models, and 
present the simulated optimal solutions. In Section 4, we validate the obtained optimal 
controls through the Hodgkin-Huxley model. 

2. Charge-Balanced Time-Optimal Control for Phase Models of Spiking 
Neuron Oscillators 

The dynamics of a periodically spiking neuron oscillator can be described by a phase 
model of the form [8] 

rlf) 

- = uj + Z(6)u(t), (1) 

where 9 denotes the phase of the oscillation, u > is neuron's natural oscillation 
frequency, and u(t) G M is the external current stimulus (control) that is applied to 
perturb the phase dynamics of the neuron. The real- valued function Z{9) is the phase 
response curve (PRC) that characterizes the infinitesimal change of the phase to an 
external control input. Conventionally, the neuron is said to spike when its phase 
6 = 2nir, where n G N. In the absence of any input u(t), the neuron spikes periodically 
at its natural frequency, while the spiking time may be advanced or delayed in a desired 
manner by the application of an appropriate weak control. 

2.1. Charge- Balanced Minimum-Time Control 

The optimal design of controls that lead to the minimum inter-spiking time of a neuron 
subject to a given bound on the control amplitude and the charge-balance constraint 
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can be formulated as a time-optimal steering problem of the form 
min T, 

u(t) 

s.t. 8 = u + Z(8)u(t), 

P = «(*), (2) 
\u(t)\ <M, Vi G [0,T], 
0(0) = 0, 0(T) = 2tt, 
p(0) = 0, p(T)=0, 

where T is the inter-spiking time required that we wish to minimize and M > is the 
bound of the control amplitude. The constraints involving the time-dependent variable 
p(t) are equivalent to the charge-balance constraint, i.e., p(t) = J^u{a)da = with 
p(0) = p(T) = 0, guaranteeing that the charge accumulated over a spiking cycle is 
zero. Note that the consideration of bounded controls is of fundamental importance 
since the phase reduction is no longer valid when the control exceeds a level that can 
be considered weak. 



2.1.1. Derivation of the Charge- Balanced Minimum- Time Control: The Hamiltonian 
of the optimal control problem as in (2) is given by 

H = X + X 1 (ou + Z(e)u) + X 2 u (3) 

where Ao > 0, Ai, and A2 are Lagrange multipliers associated with the Lagrangian, 
system dynamics, and the charge-balance constraint, respectively. According to the 
optimality conditions of the maximum principle (see Appendix A), the adjoint variables 
Ai and A 2 must satisfy the time- varying equations Ai = — 4^ and A 2 = — ^p, 
respectively, which yields 

Xx=-Xiu-^-, (4) 
A 2 = 0, (5) 

and hence A 2 is a constant. Since the Hamiltonian H is not explicitly dependent on time 
and the terminal time is free, we have H = 0, V£ G [0,T], along the optimal trajectory 
from the maximum principle. 

It is straightforward to see from a rearrangement of (3), H = Ao + X±u + (X\Z(6) + 
A 2 )w, that the control 

0< ° (6) 

minimizes the Hamiltonian H, where 

4> = X 1 Z + X 2 (7) 

is called the switching function. Hence, according to the maximum principle, u* min is a 
candidate of the optimal solution to the problem as in (2), provided = for a nonzero 
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time period does not occur. This type of controls is known as bang-bang controls, which 
takes only the extremals of the control set, e.g., — M or M in this case. The switching 
between — M and M occurs at = and the challenge is to calculate the values of 
the multipliers Ai and A2, which define the function and thus the optimal control 
sequence. 

An alternative candidate of the minimum-time control may exist. If = for some 
non-zero time interval S = [ti,t 2 ], then its derivatives 0, 0, etc., will also be equal to 
zero over S. In this case, the bang-bang control (6) may not be optimal. Such a control 
that forces the switching function and all of its derivatives to vanish over a time period 
is known as a singular control [26], and it can be calculated according to the following 
fashion. When = 0, = 0, = 0, . . ., for a given time interval S, we have 

<f ) = \ 1 Z + \ 2 = (8) 

and then, by substituting from (1), (4), and (5), the function is given by 

dZ 

* = X ^ = ( 9 ) 
which yields §§ = because u > and Ai 7^ 0. The latter is due to the non-triviality 
condition of the maximum principle, i.e., (Ao,Ai,A2) 7^ 0, since A2 = if Ai = from 
(8), which leads to Ao = from (3) as H = 0. Therefore, Ai 7^ holds along the optimal 
trajectory and §f = defines a singular trajectory, i.e., the trajectory of the system 
following a singular control. As in the calculation of 0, the second derivative can be 
obtained using (1) and |j| = to get 

<i> = x ^^( UJ + Zu )- ( 10 ) 

It is clear from (10) that if 7^ 0, the control that makes = is given by u s = —u/Z . 
In the case when = 0, we need to calculate in order to get the singular control 
u s . However, no matter how many derivatives are used, the singular control is given by 
the same form, u s = —u/Z. 

If a singular trajectory exists, then one must examine whether it is "fast" or "slow" 
compared to the bang-bang trajectory in order to determine the minimum-time control. 
Suppose that the singular control u s = —u/Z is admissible over a nonzero time interval 
S = [ti, T2]. Then, from (1) the phase velocity is equal to zero, i.e., 9 = 0, over S by the 
application of u s . This implies that the singular trajectory is slower than any feasible 
trajectory along which 9 > over S. Therefore, the charge-balanced control that spikes 
neurons in minimum time is of the bang-bang form. 

2.1.2. Computation and Synthesis of the Charge-Balanced Minimum-Time Control: 
Because the minimum spiking time of the neuron system as in (1) is achieved by a 
bang-bang control, it remains to calculate the switching points in order to synthesize 
this time-optimal control. Since = holds at the switching points, according to (8), 
these points are defined via the inverse function of the PRC, 

r-l ( ^2 



9 S = Z- L (-f i )- 
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In addition, with the Hamiltonian condition H = 0, the value of the multiplier Ai is 
given by A x = — ^ at these switching points. Without loss of generality, we let Ao = 1, 
which leads to Ai = — -. Applying this to (11) results in 

9 s = Z-\a) = Z- 1 {\ 2 u), (12) 

where A2 and u are both constants. Let Z~ l {a) have n solutions in the interval (0,27r) 
given by 81,82, ■■■ n , and define 8q = and 8 n+ i = 2ir. Then, if we start with the 
control u = M, the charge-balance constraint gives rise to the condition 

°=L*>*=xi ^j=mm de (13) 

and the total time T under this bang-bang control is represented by 

T = £/. u + (-iyz { 0)M M - (14) 



=0 



Equation (13) together with the switching conditions Z(8i) = a for % — 1, 2, . . . n define 
ri + l equations of n + 1 variables, {8\, 8 2 , ... 8 n , a}. This system of equations can be 
solved to get the set of optimal switching angles, denoted as Sm, and the constant 
a. Similarly, if we start with the control u = —M, by substituting M with — M in 
(13) we obtain the other set of solutions, denoted as S_m. The bang-bang control, 
determined by the set of switching angles, which results in the shorter spiking time is 
the charge-balanced minimum-time control, while the opposite case is a candidate for 
the charge-balanced maximum-time control. 

Alternatively, given the two sets of switching angles, the optimal switching sequence 
can be determined by computing at the switching points. We denote the vector fields 
corresponding to the constant bang controls u(t) = — M and u(t) = M by X = uj — MZ 
and Y = uj + MZ, respectively, and call the respective trajectories corresponding to 
them as X- and Y- trajectories. A concatenation of an X-trajectory followed by a Y- 
trajectory is denoted by XY , while the concatenation in the reverse order is denoted 
by YX. If < at a switching point, then the X to Y switch is optimal according to 
the switching law (6), and similarly if <fi > 0, then the Y to X switch is optimal. Since 
Ai = — l/ou at the switching points, we have 

(p = XlUJ m=-m- (15) 

Therefore, the value of ^| at the switching points defines the switching type. If ^| > 0, 
an X to Y switch is optimal and if |j| < 0, a Y to X switch is optimal. 

2.2. Charge-Balanced Maximum-Time Control 

2.2.1. (Casel: Bang-Bang Control) When the control amplitude is limited by M < 
min{|-g^y| : 8 G [0,27r)}, singular controls are not admissible since u s = —u/Z as 
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shown in Section 2.1.1. Therefore, the maximum-time control is given by the bang-bang 
form 



where <p is defined as in (7). The optimal switching sequence is determined between 
Sm and S-m, whichever results in longer spiking time. Another way to determine the 
optimal switching sequence is by evaluating ^| at the switching points as described in 
Section 2.1.2. When ^| > at a switching point, a Y to X switch is optimal, while 
when ^| < 0, an X to Y switch is optimal. 

2.2.2. (Casell: Bang-Singular-Bang Control) When singular controls are admissible, 
that is, when M > min {|^y I '■ d E [0, 2tt)}, the maximum-time control is a combination 
of bang and singular controls (see the examples in Section 3.1.2 and 3.2). The procedure 
of the optimal control synthesis is to choose a bang control that drives the system to 
a singular trajectory (a system trajectory following a singular control), staying on that 
trajectory, and then exiting at the point from which a bang control can steer the system 
to the desired terminal state. Examples involving the construction of charge-balanced 
minimum-time and maximum-time optimal controls are illustrated in Section 3. 

3. Examples 

We now apply the derived optimal control strategies to several commonly-used 
phase models characterized by various PRC's, including mathematically ideal and 
experimentally observed phase models. These examples demonstrate the applicability of 
our optimal control methods to manipulate neuron dynamics. We emphasize that these 
optimal controls are designed with respect to a given bound of the control amplitude, so 
that they can be made practical and satisfy the weak forcing assumption in the phase 
model. 

3.1. SNIPER Phase Model 

The SNIPER phase model is characterized by a type I PRC and is of the form [8] 



where u is the natural oscillation frequency of the neuron, zj, is a model-dependent 
constant, and u is the external stimulus. This model is derived from a SNIPER 
bifurcation (saddle- node bifurcation of a fixed point on a periodic orbit), which can 
be found on Type I neurons [27] such as the Hindmarsh-Rose model [28]. Neurons 
described by this model spike periodically with the natural period To = 2n/u in the 
absence of any external input u. 

Before calculating the minimum- and maximum-time spiking controls for the 
SNIPER phase model, we first examine the existence of singular trajectories. Recall 




(16) 



9 = uj + Zd (1 — cos 8)u, 



(17) 
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from (9) that the singular trajectory is defined by ^| = 0, which yields 

Zd sin 9 = 0. 

Therefore, there exist three possible singular trajectories (in this case singular points), 
9 = 0, 9 = 7r, and 9 = 2tt. The points 9 = and 9 = 2tx are infeasible singular points, 
at which the nonzero phase velocity, 9 = u, immediately forces the system away from 
these points, Hence, 9 S = ir is the only possible singular point, and the singular control 
u = —u/Z(9 s ) = —uj/(2zd) yields 9 = at 9 S , making the system stay at 9 S . 



3.1.1. Charge-Balanced Minimum-Time Control for SNIPER Phase Model: Since the 
charge-balanced minimum-time control takes the bang-bang form as shown in Section 
2.1.1, the switching points are given from (12) by 

9 K = cos- 1 (l-^\. (It 



The cosine function has two solutions in [0, 2ir) and thus there are two switching points 
9\ = 7 and 9 2 = 2ir — 7 with 7 £ [0, 7r). Because Ai = —1/u for both switching points 
and the derivative of the switching function = —Zd sin 9 < for 9 £ (0, n), if a switch 
occurs on the interval (0, n), it will be X to Y. Reversely, if a switch occurs on (n, 2tt), 
then it will be Y to X because <p > for 9 £ (tt, 2tt). It follows that an XYX trajectory 
is optimal for achieving the minimum inter-spike time. The parameter 7 that defines the 
switching points is calculated using the charge-balance constraint as in (13) by solving 
R{M,i) = 0, where 

p -m r M 

R(M,i)= / ^rjd9+ / ^rjd9. (19) 



/ u - z d (l - cos 9)M J y uj + z d (l - cos 9)M 
Then, the optimal control is given by 

f -M, 0<9<9 U 
u* mm =\ M, 9 1 <9<9 2 , (20) 
—M, 9 2 <9< 2tt, 

and by following (14) the time required to spike the neuron, namely, to reach 9 = 2tt, 
is given by 

P 4 

T = / 7^—d9. (21) 



jo u — Zd{l — cos 9)M 
Figure 1 shows the charge-balanced minimum-time control and the corresponding phase 
trajectory for the SNIPER phase model with Zd — 1, oj — 1, and M = 0.7. 

3.1.2. Charge-Balanced Maximum-Time Control for SNIPER Phase Model: There are 
two control scenarios for maximizing the spiking time of a SNIPER neuron depending 
on the control amplitude. 

(Case I: M < ^) If the bound of the control amplitude M < |^y| = | Zd(1 ^ cos g) | < 
then there exist no admissible singular controls and the maximum-time control takes 
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Figure 1. The charge-balanced minimum-time control and the corresponding phase 
trajectory for the SNIPER phase model with Zd = 1, ui = 1, and M = 0.7. 




Figure 2. The charge-balanced maximum-time control and the corresponding phase 
trajectory for the SNIPER phase model with z d = 1, ui = 1, and M = 0.4 < -f- = 0.5. 



the bang-bang form as described in Section 2.2.1. In this case, there are two switches 
and the YXY trajectory is optimal. The maximum-time control is given by 

( M, < 9 < e u 
u* max =\ —M, 9 1 <9< 9 2 , (22) 
M, 9 2 <9 < 2vr, 

where 9\ = (3, 9 2 = 2tt — (3, and the parameter (3 is obtained by solving R(—M, (3) as 
defined in (19). Figure 2 illustrates the charge-balanced maximum-time control and the 
corresponding phase trajectory for the SNIPER phase model with Zd — 1, OJ — 1, and 
M = 0.4 < ^- = 0.5. 

(Case II: M > ^-) In this case, the system can be driven along the singular trajectory 
which is optimal (slower than the bang control), and the maximum-time control takes 
the bang-singular-bang form. Because, for example, when 9 £ (0,7r), the derivative of 
the switching function <p = —z^ sin 9 < 0, and then the YX trajectory is a candidate for 
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optimality if a switch occurs. However, following an X-trajectory with u = —M < 
the singular point 9 = it is unreachable. Hence, switching in the interval (0, it) is not 
allowed, and the ^-trajectory is optimal for 9 G [0,7r). The same reasoning applies for 
the regime 9 G (tt,2tt], where F-trajectory again is optimal. As a result, the optimal 
control is of the "y-singular-y" form given by 

M, < 9 < 7T, 

6 = ^ ( 23 ) 
M, it <9 <2tt. 

Because 9 = holds along the singular trajectory (in this case the singular point 9 S = it), 
the time duration over which the system stays on it is calculated according to the charge- 
balance constraint. Let t\ and t2 denote the times for which the first bang control and 
the singular control are applied, respectively. Since ti is the time that the system takes 
to reach 9 S = it by a V-trajectory, we have 

h = / 7 7TT—d9. (24) 

By symmetry, the amount of time that the system takes following a ^-trajectory from 
9 = it to 9 = 2tt is also t±. Then, ti is given by 

_ AMzah 

£2 — 

00 

in order to fulfill the charge-balance constraint. Now the charge-balanced maximum- 
time control can be stated in terms of time as 

M, < t < ti, 

-4, ti<t<ti + t 2 , (25) 
M, t 1 + t 2 <t<t 2 + 2t 1 . 

Figure 3 shows the maximum-time charge-balanced control and the corresponding phase 
trajectory for the SNIPER phase model with Zd = 1, to = 1, and M = 0.7 > = 0.5. 

In the following, we demonstrate the robustness of our analytical method to 
construct optimal controls for spiking neurons of arbitrary practical PRCs through the 
Hodgkin-Huxley and Morris-Lecar phase models. 



3.2. Hodgkin-Huxley Phase Model 

The Hodgkin-Huxley model is a nonlinear system that characterizes the propagation 
and initiation of the action potential in a squid axon [23]. For the set of parameter 
values given in [8], the system exhibits periodic motion with natural frequency u = 
0A3rad/ms. Its PRC and the first and second derivatives of the PRC are depicted 
in Figure 4(a). To proceed the calculations, we approximate the numerically obtained 
PRC with eight harmonic terms given by 

8 

Z(0) = J]a < sin(6 i + c i ), (26) 

i=l 
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Figure 3. The maximum-time charge-balanced control and corresponding phase 
trajectory for the SNIPER phase model with Zd = 1,oj = 1, and M = 0.7. 
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Table 1. The coefficients of the equation (26) for the Hodgkin-Huxley PRC. 



where the coefficients a,j, bi and c, are obtained by least squares fit and given in Table 1. 
In this case, there are two possible singular points, 6 St x = 3.34 and 9 S ^ = 4.58, satisfying 
dZ(6)/d6 = 0. 

The charge-balanced minimum-time control, which is of the YXY form, and the 
resulting phase trajectory for the control amplitude bound M = 0.7fiAcm~ 2 are shown 
in Figure 4(b). The charge-balanced maximum-time controls can take the bang-bang or 
the bang-sigular-bang form depending on the bound M. The cases for M = 0.7 ixAcvrc' 2 
and M = 3.0 fiAcm~ 2 are illustrated in Figure 4(c) and Figure 4(d), respectively. The 
detailed derivations of these optimal controls are presented in Appendix B. 

3.3. Morris-Lecar Phase Model 

The Morris-Lecar neuron model was originally developed to capture the oscillatory 
behavior of barnacle muscle fibers [25]. It has been observed through experiments that 
the PRC for an Aplysia motoneuron is extremely similar to that of a Morris-Lecar PRC 
[29]. We consider a Morris-Lecar system with parameter values given in [12], which has 
a natural frequency u = 0.283rad/ms. The PRC is approximated by (26) with the 
coefficients shown in Table 2 and is illustrated, with its derivatives, in Figure 5(a). 

Three examples are made to show the different structures of the optimal controls 
that are associated with various values of M for the Morris-Lecar phase model. The 
charge-balanced minimum-time control and the resulting phase trajectory for M = 
0.01 fiAcm~ 2 are given in Figure 5(b). The charge-balanced maximum-time controls 
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Figure 4. 4(a) The Ho dgkin- Huxley PRC Z{6) and its derivatives, |f and ^gf . 4(b) 
The charge-balanced minimum-time control and the corresponding phase trajectory for 
the Hodgkin-Huxley phase model with respect to the bound on the control amplitude 
M = 0.7 /j,Acm~ 2 . 4(c) and 4(d) show the charge-balanced maximum-time controls 
and the corresponding phase trajectories for M = 0.7 nAcm~ 2 and M = 3.0 \iAcmT 2 . 
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Table 2. The coefficients of the equation (26) for the Morris-Lecar PRC 



and the respective trajectories subject to M — 0.005 fiAcm~ 2 and M = 0.04 fiAcm~ 2 
are given in Figure 5(c) and Figure 5(d), respectively. The derivations of these optimal 
controls follow a similar procedure presented in Appendix B. 

4. Validation of Phase Model Reduction to Full State-Space Model 

Because phase models of importance to applications are reductions of original higher 
dimensional state-space systems, we explore in this section the extent to which controls 
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Figure 5. 5(a) The Morris-Lccar PRC Z(0) and its derivatives, ^§ and 5(b) The 
charge-balanced minimum-time control and the corresponding phase trajectory for the 
Morris-Lecar phase model with respect to the bound on the control amplitude M = 
0.01 /j,Acm~ 2 . 5(c) and 5(d) show the charge-balanced maximum-time controls and 
the corresponding phase trajectories with M = 0.005 /j,Acm~ 2 and M = 0.04 fj,Acm~ 2 , 
respectively. 



synthesized using the former can achieve the desired objectives when applied to the 
latter. This will provide insight into the limits of the model reduction with respect to 
control synthesis, and allow the relationship to be calibrated for practical applications 
where the weak forcing assumption must be relaxed. Such an important validation is 
largely lacking in the literature. 

We validate our optimal control strategies derived based on the phase models with 
the corresponding original state-space models. Specifically, we consider the Hodgkin- 
Huxley model. Note that an analytical derivation of the optimal controls directly 
from the state-space models is in general intractable and computationally expensive. 
A validation of the minimum and maximum spiking times with respect to the bound 
on the control amplitude is depicted in Figure 6, where the feasible spiking times are 
indicated as the shaded area. Each asterisk point on this graph represents the Hodgkin- 
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Figure 6. A characterization of the realizable spiking times with respect to the 
bound on the control amplitude, M € [0,2.5], for the Hodgkin-Huxley phase model. 
The shaded region indicates the feasible spiking range resulting from the minimum- 
and maximum-time controls. Those minimum times (left to the natural spiking time 
To = 14.6 ms) are obtained by YXY controls and maximum-times (right to To) are 
obtained by XYX, YXY and Y-singular-F controls depending on M . 

Huxley neuron spiking time achieved by a particular form of the optimal control. The 
points correspond to minimum spiking times, which are less than the natural spiking 
time To = 14.64 ms, are obtained by YXY controls, whereas the points correspond to 
maximum spiking times may be obtained by three structurally different controls, i.e., 
XYX, YXY, and F-singular-y controls. This figure illustrates the limits on possible 
spiking times of the Hodgkin-Huxley model, which is important to the design of practical 
control inputs. For example, the knowledge of the feasible spiking range is helpful in 
designing optimal controls with other objectives such as minimum power controls [13]. 

The optimal controls derived based on the Hodgkin-Huxley phase model, shown in 
Figure 4(b) and 4(c), are applied to the full Hodgkin-Huxley model, and a repeated 
application of such controls results in the desired spiking trains as displayed in 
Figure 7(a) and 7(b). The respective minimum and maximum spiking times induced 
from these optimal controls subject to the amplitude bound M = 0.7 fiAcm~ 2 are 
13.5 ms and 16.37 ms in the phase model and 13.65 ms and 17.13 ms in the full state- 
space model. Such an inconsistence is due to the model reduction, however, the resulting 
spiking behavior of the full Hodgkin-Huxley model shows great qualitative agreement 
with that of the phase model. The variation of the absolute errors between the actual 
and designed spiking times is shown in Figure 8, where the spiking behavior predicted 
based on the phase model matches better the full state-space model towards the weak 
forcing region. 
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(a) Uncontrolled and controlled spiking trains for minimum time with amplitude 
M = 0.7 \xAcmT 2 of Hodgkin-Huxley neuron. 
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(b) Uncontrolled and controlled spiking trains for maximum time with amplitude 
M = 0.7 fiAcm~ 2 of Hodgkin-Huxley neuron. 

Figure 7. Application of derived optimal controls according to phase models to full 
Hodgkin-Huxley model 




Spiking Time According to Phase Model (ms) 



Figure 8. The absolute error in the spiking time when applying the charge-balanced 
time-optimal controls derived based on the Hodgkin-Huxley phase model to its full 
state-space model. The bound of the control amplitude is indicated as the color bar. 



5. Conclusion and Future Work 

In this paper, we investigated time-optimal controls for phase models of spiking neuron 
oscillators. In particular, we derived charge-balanced controls that lead to the minimum 
and the maximum inter-spike time of a neuron for a given bound on the control 
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amplitude. We showed that such optimal controls involve bang-bang and bang-singular- 
bang structures depending on the allowable control amplitude. Although the amplitude 
level of weak forcing in the phase model is not practically quantifiable and can be greatly 
dependent on the dynamics of the system, our optimal control solutions were constructed 
for an arbitrary choice of bounds on the control amplitude, which accounts for this 
practical issue. We apply the derived optimal spike timing controls to commonly-used 
phase models of oscillatory neurons to demonstrate their applicability to neuroscience. 
The methodology presented in this paper is general and can be applied not only to 
oscillatory neuron systems, but also to any oscillating system that can be represented 
by phase models including biological, chemical, electrical, and mechanical oscillators. 

The theoretical results of this work characterize the fundamental limits on neuron 
spiking times that can be achieved by use of a charge-balanced bounded external 
input, and have potential impact on the improvement and development of innovative 
therapeutic procedures for neurological disorders. The extension of this work to the 
optimal control of a neuron population is of fundamental and practical importance. 
Our recent work has shown that an ensemble of uncoupled neurons is controllable and 
the minimum-power controls that spike a network of heterogeneous neurons can be 
constructed using a multidimensional pseudospectral method [19]. We plan to extend 
this current work to investigate the controllability and optimal controls of a network of 
coupled neurons. 



Appendix A. The Pontryagin's Maximum Principle 

Theorem 1 (Time-Optimal Control [24]) Let u*{t)) be a time-optimal 

controlled trajectory that transfers the initial condition x(0) = xq into the terminal 
state x(T) = xt ■ Then, it is a necessary condition for optimality that there exists a 
constant Ao > and nonzero, absolutely continuous row vector function X(t) such that: 

(i) A satisfies the so-called adjoint equation 

OH 

A(t) = -—(\ Q ,\(t),x*(t),u*(t)) 

(ii) For < t < T the function u y H(\q, X(t), u) attains its minimum over the 
control set U at u = u*(t). 

(hi) H{X ,X(t),x*{t),u*(t)) = 0. 



Appendix B. The Derivation of Time-Optimal Controls for the 
Hodgkin-Huxley Phase Model 

Charge-Balanced Minimum-Time Control for Hodgkin-Huxley Phase Model 

The Hodgkin-Huxley PRC given in Figure 4(a) has at most two singular trajectories 
(points), 6> Sj i = 3.34 and 9 Si2 = 4.58, calculated by the condition 9z g ^ = 0. According 
to the shape of this PRC, there exist at most two switching points satisfying Z{9) = a, 
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where a is a constant defined in (12). Since the minimum-time control takes the bang- 
bang form as shown in Section 2.1.1, it requires to calculate the switching points and 
determine the type of the switching at these points for the optimal control synthesis. 
At the switching points, <fi = —dZ/89 is given by (15), and hence a Y to X switch 
may occur in the region R\ = [0,8 s .i] or R% = [6 s ^,2tt], and an X to Y switch may 
occur in R 2 = [0 s ,i, #3,2] • This implies that bang-bang controls with one switch, such as 
the XY or YX form, are not feasible solutions because these controls will violate the 
charge-balance constraint. Consequently, the optimal control has two switching points, 
and the candidate is either a YXY trajectory with one switch in the interval R\ and 
one in R2, or an XYX trajectory with one switch in R2 and one in R3. We can further 
simplify the possible intervals of switching by observing the shape of the PRC. The 
Hodgkin-Huxley PRC depicted in Figure 4(a) has three zeros at 9 r> i = 0, 9 r> 2 = 3.86, 
and 9 r $ = 2ir. Therefore, for an optimal YXY trajectory the first and the second switch 
will occur in [0, 8 St i] and [^,1,^,2], respectively, and for an optimal XYX trajectory, 
they will occur in [0 r ,2, ^,2] and [# S) 2,27r], respectively. The minimum-time control is 
then selected between these two. Note that for a given bound M, it may not be possible 
to have both XYX and YXY solutions. For example, if the bound is M = 0.7, then 
the only feasible optimal solution is YXY. In this case, the two switching points Q\ and 
62 can be calculated through 



= J ^TWW) de + L u-Mzyf + k ^TMW) d9ABA) 

Z(9 1 ) = Z(9 2 ), (B.2) 

and the control is then given by 

M, < 6 < 0x, 

—M, e 1 < e < 6 2 , 

M, 6 2 < 9 < 2tt. 

Charge- Balanced Maximum- Time Control for Hodgkin-Huxley Phase Model 

In the case of the maximum-time control, the two singular points, 6 Si \ and 9 Sj 2, are 
candidates for the optimal trajectory because they are slower than the bang trajectories 
as proved in Section 3.1.2. Letting 9 = in (1), we find the controls that keep the 
trajectory at the singular points are u s \ — — >, = 3.50 and u s2 = — , = —2.15. 
There exist three cases when constructing maximum-time controls according to M and 
thus to the feasibility of u s> \ and u Sj2 . 

( Case I: M < \u s>2 \) In this case, both the singular points Sj i and 9 Sy2 are infeasible. 
Therefore, the optimal control is bang-bang and is in fact the opposite of the minimum- 
time control described above. Similar to the minimum-time case, we can calculate 
the corresponding XYX and YXY solutions and choose the maximum time achieved 
between these scenarios. For example, consider the bound M = 0.7, then the only 
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solution is XYX and the two switching points are calculated by substituting M with 
— M in (B.l) and solving (B.l) and (B.2). The optimal bang-bang control is then given 

by 

' —M, < 9 < X , 

m, e 1 < e < e 2 , 

-M, 9 2 <9 < 2vr. 

( Case II: |w S)2 | < M < \u s< i\) In this case, 9 Sj2 is the only feasible singular trajectory 
(point) generated by the singular control u S)2 = —2.15 < 0. Because there are only two 
switching points allowed in the optimal trajectory, this together with the fact that u Sj2 
is of negative charge forces the optimal control to take the "y-singular-y" form given 
by 

o < 9 < e 8t2 , 

9 = #s,2, 

9 S ,2 < 9 < 2vr. 

Similar to the SNIPER phase model described in Section 3.1.2, the time it takes to 
reach the singular point is given by, 

tl = Jo u + Z(9)M d9 
and the time required to reach the target point 2n from the point 9 S>2 is 

h = L u 1 + Z(6)M d °- 

The time during which the trajectory stays on 9 s 2 is determined by the charge-balance 
constraint and is given by 

{h + h)M 




h 

U s ,2 

Now, the optimal control can be stated in terms of time as 

M, < t < ti, 

u s ,2, h<t< h+t 2 , (B.3) 
M, tx + t 2 <t<tx + t 2 + 1 3 . 

(Case III: \u Sj i\ < M) In this case, staying on either singular point is possible by 
using an appropriate control. Furthermore, since the two singular controls have opposite 
signs, the charge-balance constraint can be preserved by staying for an appropriate 
time period at each singular point. As a result, the spiking time can be arbitrarily 
delayed, which may not be of practical interest due to the requirement of relatively high 
amplitude. 
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